Libraries:

# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)
library(ggeffects)
library(performance)

Set working directory

knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')

Source sapling data:

source("Scripts/SA_Import.R")

Pivot wider to create dataframe where each row is for one plot and has total details for each species group

sa_merge2 <- sapling_merge %>% 
  pivot_wider(names_from = Species_Groups, values_from = Total_Tally)

Import time since data and add it to the PIRI dataset

time_since <- read_csv("CleanData/Treat_Year_Data.csv")

sa_merge3 <- merge(sa_merge2, time_since, by = 'Site')
#log transform time from treatment data
sa_merge3$l.TFT <- log(sa_merge3$Time_from_Treat)

Run the ‘Add_BA’ script and merge with dataset:

source("Scripts/Add_BA.R")

# merge with sapling dataset -------------------
sa_merge4 <- merge(sa_merge3, prism_BA, by = 'Plot_No')

Run ‘Ground_Data.R’ script and add it to sapling dataset:

source("Scripts/Ground_Data.R")

# merge with sapling dataset -------------------
sa_merge5 <- merge(sa_merge4, ground3, by = 'Plot_No')

Source and add veg data:

source("Scripts/Veg_Data.R")

# merge with sapling dataset -------------------
sa.all <- merge(sa_merge5, veg3, by = "Plot_No")

rm(sa_merge3,
   sa_merge4,
   sa_merge5)

Sapling count data is taken in 25m2 plots; basal area is measured in hectares; veg and soil data is taken a 1m2 plots.

Data will be standardized for comparisons at multiple scales

Create a dataframe keeping sapling at 25m2 observation levels & log transforming them

sa.all3 <- sa.all

sa.all3$piri.ba_s <- scale(sa.all3$PIRI.BA_HA)
sa.all3$avgld_s <- scale(sa.all3$avgLD)
sa.all3$vegt_s <- scale(sa.all3$Veg_Total)
sa.all3$min_s <- scale(sa.all3$Mineral_Soil)
sa.all3$so_s <- scale(sa.all3$Shrub_Oak)
sa.all3$other_s <- scale(sa.all3$Other)
sa.all3$ba.ha_s <- scale(sa.all3$BA_HA)
sa.all3$piri_s <- scale(sa.all3$PIRI)

sa.all3 <- sa.all3 %>% 
  select(Treat_Type, Region, Site, Plot_No, PIRI, piri_s, Shrub_Oak, so_s, Other, other_s, Time_from_Treat, l.TFT, BA_HA, ba.ha_s, PIRI.BA_HA, piri.ba_s, Mineral_Soil, min_s, Litter_Duff, avgLD, avgld_s, Veg_Total, vegt_s) %>% 
  arrange(Treat_Type)

Looking at pairwise with sa.all3 dataset (sapling counts at 25m2 observation level)

#not log transformed (25m2)
sa3.num <- sa.all3 %>% 
  select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)

ggpairs(sa3.num)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggpairs(sa3.num, aes(color = Treat_Type))
## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

rm(sa3.num)

No strong correlations

major personal breakthrough for dummies. For zero inflation, lets see if any TT or region is lacking alltogether in PIRI or SO saplings

sa.all4 <- sa.all3 %>% 
  group_by(Region)

tapply(sa.all4$PIRI, sa.all4$Region, summary)
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.02459 0.00000 2.00000 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2013  0.0000  9.0000 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.07273 0.00000 1.00000 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.09574 0.00000 3.00000
# no PIRI in LI,
tapply(sa.all4$PIRI, sa.all4$Treat_Type, summary) 
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2564  0.0000  9.0000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.04902 0.00000 2.00000 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.05357 0.00000 2.00000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.05882 0.00000 2.00000 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01149 0.00000 1.00000

Now that I know that there are no PIRI in LI, I want to start model elimination again using ziformula = ~ Region

I have previously done model elimination with zero inflated poisson distro & with model w.o TT. Neither are better than zi nbinom2 models. zi ~1 has distribution issues, ~Region does not

Model has much better fit and zero inflation issues are resolved when treatment type variable is included. Time to do elimination again but with TT; all models fail with nbinom1 distro, but works with nbinom2 (these notes are from when I was using a zi poisson distro, using nbinom TT no longer seems to be important to model fit, but we will keep for the overarching story)

Checking fixed effects for collinearity and looking at general summary stats to make sure they aren’t 0/1 or same/function of other category

summary(sa.all3)
##   Treat_Type           Region              Site              Plot_No      
##  Length:498         Length:498         Length:498         Min.   :  70.0  
##  Class :character   Class :character   Class :character   1st Qu.: 359.5  
##  Mode  :character   Mode  :character   Mode  :character   Median : 614.0  
##                                                           Mean   : 613.5  
##                                                           3rd Qu.: 887.5  
##                                                           Max.   :1143.0  
##       PIRI              piri_s.V1        Shrub_Oak            so_s.V1      
##  Min.   :0.00000   Min.   :-0.157452   Min.   :0.0000   Min.   :-0.273816  
##  1st Qu.:0.00000   1st Qu.:-0.157452   1st Qu.:0.0000   1st Qu.:-0.273816  
##  Median :0.00000   Median :-0.157452   Median :0.0000   Median :-0.273816  
##  Mean   :0.09438   Mean   : 0.000000   Mean   :0.2691   Mean   : 0.000000  
##  3rd Qu.:0.00000   3rd Qu.:-0.157452   3rd Qu.:0.0000   3rd Qu.:-0.273816  
##  Max.   :9.00000   Max.   :14.857482   Max.   :7.0000   Max.   : 6.849481  
##      Other            other_s.V1      Time_from_Treat     l.TFT      
##  Min.   :0.0000   Min.   :-0.345712   Min.   : 1.0    Min.   :0.000  
##  1st Qu.:0.0000   1st Qu.:-0.345712   1st Qu.: 3.0    1st Qu.:1.099  
##  Median :0.0000   Median :-0.345712   Median : 5.0    Median :1.609  
##  Mean   :0.3514   Mean   : 0.000000   Mean   :10.9    Mean   :1.804  
##  3rd Qu.:0.0000   3rd Qu.:-0.345712   3rd Qu.: 8.0    3rd Qu.:2.079  
##  Max.   :9.0000   Max.   : 8.508476   Max.   :33.0    Max.   :3.497  
##      BA_HA            ba.ha_s.V1        PIRI.BA_HA      piri.ba_s.V1    
##  Min.   : 0.00   Min.   :-1.4237625   Min.   : 0.0   Min.   :-1.285691  
##  1st Qu.: 7.00   1st Qu.:-0.8185636   1st Qu.: 5.0   1st Qu.:-0.775675  
##  Median :16.00   Median :-0.0404508   Median :11.0   Median :-0.163656  
##  Mean   :16.47   Mean   : 0.0000000   Mean   :12.6   Mean   : 0.000000  
##  3rd Qu.:23.00   3rd Qu.: 0.5647481   3rd Qu.:18.0   3rd Qu.: 0.550367  
##  Max.   :51.00   Max.   : 2.9855436   Max.   :44.0   Max.   : 3.202451  
##   Mineral_Soil          min_s.V1        Litter_Duff        avgLD       
##  Min.   : 0.0000   Min.   :-0.197640   Min.   : 0.00   Min.   : 0.000  
##  1st Qu.: 0.0000   1st Qu.:-0.197640   1st Qu.:90.50   1st Qu.: 2.250  
##  Median : 0.0000   Median :-0.197640   Median :90.50   Median : 3.875  
##  Mean   : 0.7199   Mean   : 0.000000   Mean   :88.78   Mean   : 4.179  
##  3rd Qu.: 0.0000   3rd Qu.:-0.197640   3rd Qu.:90.50   3rd Qu.: 5.500  
##  Max.   :50.5000   Max.   :13.666913   Max.   :90.50   Max.   :13.250  
##      avgld_s.V1        Veg_Total           vegt_s.V1     
##  Min.   :-1.652607   Min.   :  4.00   Min.   :-1.722482  
##  1st Qu.:-0.762793   1st Qu.: 37.50   1st Qu.:-0.738485  
##  Median :-0.120151   Median : 58.00   Median :-0.136337  
##  Mean   : 0.000000   Mean   : 62.64   Mean   : 0.000000  
##  3rd Qu.: 0.522492   3rd Qu.: 83.38   3rd Qu.: 0.609004  
##  Max.   : 3.587405   Max.   :202.00   Max.   : 4.093383
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
findLinearCombos(sa.all3[,4:23]) #shows no linear combos. Don't get how the rank deficiency is happening
## $linearCombos
## $linearCombos[[1]]
## [1] 5 2 3 4
## 
## $linearCombos[[2]]
## [1] 7 2 3 6
## 
## $linearCombos[[3]]
## [1] 11  2  3 10
## 
## $linearCombos[[4]]
## [1] 13  2  3 12
## 
## $linearCombos[[5]]
## [1] 15  2  3 14
## 
## $linearCombos[[6]]
## [1] 18  2  3 17
## 
## $linearCombos[[7]]
## [1] 20  2  3 19
## 
## 
## $remove
## [1]  5  7 11 13 15 18 20
tapply(sa.all3$Other, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.3443  0.0000  6.0000 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   0.411   0.000   6.000 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2143  0.0000  6.0000 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.6727  0.0000  9.0000 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.3511  0.0000  4.0000
tapply(sa.all3$Shrub_Oak, sa.all3$Region, summary) #none in LI
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## $MA
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.000000 0.000000 0.006494 0.000000 1.000000 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   1.964   3.000   7.000 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2447  0.0000  6.0000
tapply(sa.all3$BA_HA, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   11.00   14.98   23.00   51.00 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   11.00   21.00   19.81   28.00   39.00 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    9.00   14.00   15.32   21.00   44.00 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   13.50   23.00   22.13   31.00   44.00 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   11.00   14.37   21.00   48.00
tapply(sa.all3$PIRI.BA_HA, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       2       9      10      16      41 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   11.00   16.00   17.41   23.00   37.00 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00    9.00   11.44   16.00   34.00 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    8.00   18.00   19.18   30.00   44.00 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00    9.00   10.32   16.00   32.00
tapply(sa.all3$Veg_Total, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.00   51.50   71.00   70.96   91.88  166.50 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    7.00   23.00   39.00   44.52   61.00  128.00 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    6.50   33.75   52.75   56.93   74.38  161.50 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   14.50   41.25   64.00   62.97   79.75  114.00 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.00   46.00   62.25   75.09   96.62  202.00
tapply(sa.all3$Mineral_Soil, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   1.586   0.500  30.500 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   0.589   0.500  15.500 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.6623  0.0000 50.5000 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01818 0.00000 0.50000 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2021  0.0000  8.0000
tapply(sa.all3$avgLD, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.000   1.750   1.977   2.750   7.250 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.75    4.00    5.25    5.61    6.75   12.50 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.250   3.562   5.000   5.347   6.750  13.250 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.750   2.500   3.500   3.686   4.500   9.750 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.750   2.812   4.000   4.298   5.438  12.250
tapply(sa.all3$Other, sa.all3$Treat_Type, summary) #all present
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.8889  1.0000  6.0000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.4118  0.0000  9.0000 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.1429  0.0000  2.0000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.06618 0.00000 2.00000 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.1379  0.0000  2.0000
tapply(sa.all3$Shrub_Oak, sa.all3$Treat_Type, summary) #none in mowrx or springrx
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.7094  0.0000  7.0000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.4902  0.0000  6.0000 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01786 0.00000 1.00000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0
tapply(sa.all3$BA_HA, sa.all3$Treat_Type, summary) #all present
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   18.00   28.00   26.78   37.00   51.00 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    9.00   16.00   17.44   23.00   44.00 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   11.00   12.11   18.00   34.00 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   7.000   9.059  14.000  37.000 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   10.00   14.00   15.85   21.00   41.00
tapply(sa.all3$PIRI.BA_HA, sa.all3$Treat_Type, summary) #all present
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    9.00   18.00   17.05   25.00   41.00 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   11.00   14.45   20.25   44.00 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   5.000   9.000   9.804  14.500  30.000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   7.000   7.816  11.000  37.000 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    9.00   14.00   13.75   18.00   34.00
tapply(sa.all3$Veg_Total, sa.all3$Treat_Type, summary) #all present
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    7.00   26.50   46.00   46.85   61.00  115.50 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    14.5    46.0    68.0    72.7    91.0   202.0 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.00   36.12   55.00   55.21   68.62  132.00 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.00   37.00   66.75   70.47   94.38  196.50 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   16.00   49.50   61.50   64.65   81.00  131.00
tapply(sa.all3$Mineral_Soil, sa.all3$Treat_Type, summary) #all present
## $Control
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.000000 0.000000 0.008547 0.000000 0.500000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   0.201   0.000   8.000 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   1.652   0.000  50.500 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   1.199   0.500  30.500 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.9368  0.5000 15.5000
tapply(sa.all3$avgLD, sa.all3$Treat_Type, summary) #all present
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.500   3.250   4.750   5.024   6.750  12.500 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.750   3.250   4.375   4.529   5.750   9.750 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.250   3.688   4.775   5.291   6.750  11.000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.500   2.750   2.816   3.812  13.250 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.500   1.500   3.500   4.046   5.250  11.500

PIRI Model

Best Model AIC: 245.8

Running model elimination used standardized effects and zi by Region

sa.m1s <- glmmTMB(PIRI ~ Treat_Type + piri.ba_s + avgld_s + vegt_s + min_s + so_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; false convergence (8). See vignette('troubleshooting'),
## help('diagnose')
#won't converge

#test piri ba
sa.m2s <- glmmTMB(PIRI ~ Treat_Type  + avgld_s + vegt_s + min_s + so_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m2s) #248
## [1] 247.957
# test ba/ha
sa.m3s <- glmmTMB(PIRI ~ Treat_Type  + avgld_s + vegt_s + min_s + so_s + other_s  + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m3s) #251.7
## [1] 251.6733
lrtest(sa.m2s, sa.m3s) #keep ba
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s + 
##     ba.ha_s + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s + 
##     offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1  19 -104.98                       
## 2  18 -107.84 -1 5.7163    0.01681 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test SO
sa.m4s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s  + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m4s) #248
## [1] 248.2719
lrtest(sa.m4s, sa.m2s) #drop
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s + 
##     offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s + 
##     ba.ha_s + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  18 -106.14                     
## 2  19 -104.98  1 2.3148     0.1281
# test other
sa.m5s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m5s) #251.1
## [1] 251.0924
lrtest(sa.m4s, sa.m5s) #keep other
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s + 
##     offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + offset(l.TFT) + 
##     (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1  18 -106.14                       
## 2  17 -108.55 -1 4.8206    0.02812 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test mineral
sa.m6s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + other_s  + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m6s) #247
## [1] 247.0116
lrtest(sa.m4s, sa.m6s) #drop
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s + 
##     offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + other_s + offset(l.TFT) + 
##     (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  18 -106.14                     
## 2  17 -106.51 -1 0.7397     0.3897
# test veg
sa.m7s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s  + other_s  + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m7s) #245.8
## [1] 245.8436
# test ld
sa.m8s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + other_s  + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m8s) #244.8 
## [1] 244.7564
sa.m8s_sr <- simulateResiduals(sa.m8s, n = 1000, plot = TRUE) #passes

testResiduals(sa.m8s_sr)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.029239, p-value = 0.7881
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.4242, p-value = 0.826
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00072289156626506 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.029239, p-value = 0.7881
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.4242, p-value = 0.826
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00072289156626506 ) 
##                                                  0
# passes


emmeans(sa.m8s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response') # the model is rank deficient
## $emmeans
##  Treat_Type response  SE  df asymp.LCL asymp.UCL
##  Control      0.0272 NaN Inf       NaN       NaN
##  FallRx       0.0235 NaN Inf       NaN       NaN
##  Harvest      0.1676 NaN Inf       NaN       NaN
##  MowRx        0.1509 NaN Inf       NaN       NaN
##  SpringRx     0.0193 NaN Inf       NaN       NaN
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio  SE  df null z.ratio p.value
##  Control / FallRx   1.157 NaN Inf    1     NaN     NaN
##  Control / Harvest  0.162 NaN Inf    1     NaN     NaN
##  Control / MowRx    0.180 NaN Inf    1     NaN     NaN
##  Control / SpringRx 1.409 NaN Inf    1     NaN     NaN
##  FallRx / Harvest   0.140 NaN Inf    1     NaN     NaN
##  FallRx / MowRx     0.156 NaN Inf    1     NaN     NaN
##  FallRx / SpringRx  1.217 NaN Inf    1     NaN     NaN
##  Harvest / MowRx    1.111 NaN Inf    1     NaN     NaN
##  Harvest / SpringRx 8.678 NaN Inf    1     NaN     NaN
##  MowRx / SpringRx   7.812 NaN Inf    1     NaN     NaN
## 
## P value adjustment: tukey method for comparing a family of 2 estimates 
## Tests are performed on the log scale
sa.m9s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m9s) #245.9 with veg; 245.8 with ld
## [1] 245.8436
summary(sa.m9s)
##  Family: nbinom2  ( log )
## Formula:          
## PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) +  
##     (1 | Site/Plot_No)
## Zero inflation:        ~Region
## Data: sa.all3
## 
##      AIC      BIC   logLik deviance df.resid 
##    245.8    313.2   -106.9    213.8      482 
## 
## Random effects:
## 
## Conditional model:
##  Groups       Name        Variance      Std.Dev. 
##  Plot_No:Site (Intercept) 0.00000003136 0.0001771
##  Site         (Intercept) 1.00553141354 1.0027619
## Number of obs: 498, groups:  Plot_No:Site, 498; Site, 47
## 
## Dispersion parameter for nbinom2 family (): 0.332 
## 
## Conditional model:
##                    Estimate Std. Error z value       Pr(>|z|)    
## (Intercept)         -5.4223     0.8811  -6.154 0.000000000755 ***
## Treat_TypeFallRx    -0.2051     0.9401  -0.218         0.8273    
## Treat_TypeHarvest    1.8988     1.3009   1.460         0.1444    
## Treat_TypeMowRx      1.6541     1.0395   1.591         0.1116    
## Treat_TypeSpringRx  -0.2551     1.4500  -0.176         0.8603    
## ba.ha_s              0.8721     0.3623   2.407         0.0161 *  
## other_s             -1.1511     0.6811  -1.690         0.0910 .  
## avgld_s             -0.2776     0.2909  -0.954         0.3398    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)     1.521      1.019   1.492    0.136
## RegionLI       22.406  53458.279   0.000    1.000
## RegionMA      -20.529  22550.145  -0.001    0.999
## RegionME       -2.795      3.323  -0.841    0.400
## RegionNH       -1.687      1.457  -1.158    0.247

Best model (lowest AIC without rank deficiency) appears to be model sa.m6 with TT, piri BA, agvLD, Veg total. Test model fit:

Model without avgld is rank deficient

sa.m9s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m9s) #245.9 with veg; 245.8 with ld
## [1] 245.8436
summary(sa.m9s)
##  Family: nbinom2  ( log )
## Formula:          
## PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) +  
##     (1 | Site/Plot_No)
## Zero inflation:        ~Region
## Data: sa.all3
## 
##      AIC      BIC   logLik deviance df.resid 
##    245.8    313.2   -106.9    213.8      482 
## 
## Random effects:
## 
## Conditional model:
##  Groups       Name        Variance      Std.Dev. 
##  Plot_No:Site (Intercept) 0.00000003136 0.0001771
##  Site         (Intercept) 1.00553141354 1.0027619
## Number of obs: 498, groups:  Plot_No:Site, 498; Site, 47
## 
## Dispersion parameter for nbinom2 family (): 0.332 
## 
## Conditional model:
##                    Estimate Std. Error z value       Pr(>|z|)    
## (Intercept)         -5.4223     0.8811  -6.154 0.000000000755 ***
## Treat_TypeFallRx    -0.2051     0.9401  -0.218         0.8273    
## Treat_TypeHarvest    1.8988     1.3009   1.460         0.1444    
## Treat_TypeMowRx      1.6541     1.0395   1.591         0.1116    
## Treat_TypeSpringRx  -0.2551     1.4500  -0.176         0.8603    
## ba.ha_s              0.8721     0.3623   2.407         0.0161 *  
## other_s             -1.1511     0.6811  -1.690         0.0910 .  
## avgld_s             -0.2776     0.2909  -0.954         0.3398    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)     1.521      1.019   1.492    0.136
## RegionLI       22.406  53458.279   0.000    1.000
## RegionMA      -20.529  22550.145  -0.001    0.999
## RegionME       -2.795      3.323  -0.841    0.400
## RegionNH       -1.687      1.457  -1.158    0.247
sa.m9s_sr <- simulateResiduals(sa.m9s, n = 1000, plot = TRUE) #passes

testResiduals(sa.m9s_sr)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.030633, p-value = 0.7383
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.30212, p-value = 0.862
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000883534136546185 ) 
##                                                   0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.030633, p-value = 0.7383
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.30212, p-value = 0.862
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000883534136546185 ) 
##                                                   0
# passes

Pairwise

emmeans(sa.m9s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
##  Treat_Type response     SE  df asymp.LCL asymp.UCL
##  Control      0.0268 0.0236 Inf   0.00477     0.151
##  FallRx       0.0219 0.0206 Inf   0.00345     0.139
##  Harvest      0.1792 0.1821 Inf   0.02444     1.314
##  MowRx        0.1403 0.1040 Inf   0.03279     0.600
##  SpringRx     0.0208 0.0280 Inf   0.00148     0.292
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio     SE  df null z.ratio p.value
##  Control / FallRx   1.228  1.154 Inf    1   0.218  0.9995
##  Control / Harvest  0.150  0.195 Inf    1  -1.460  0.5889
##  Control / MowRx    0.191  0.199 Inf    1  -1.591  0.5030
##  Control / SpringRx 1.291  1.871 Inf    1   0.176  0.9998
##  FallRx / Harvest   0.122  0.166 Inf    1  -1.542  0.5352
##  FallRx / MowRx     0.156  0.170 Inf    1  -1.704  0.4314
##  FallRx / SpringRx  1.051  1.577 Inf    1   0.033  1.0000
##  Harvest / MowRx    1.277  1.562 Inf    1   0.200  0.9996
##  Harvest / SpringRx 8.619 14.264 Inf    1   1.302  0.6903
##  MowRx / SpringRx   6.748 10.028 Inf    1   1.285  0.7008
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log scale

None significantly different from eachother

Plot model

sa.piri1 <- ggpredict(sa.m9s, terms = c("avgld_s", "Treat_Type"))
#control and fallrx are almost completely identical - if I wanted to nudge/jitter
# sa.piri2 <- sa.piri1 %>% 
#   filter(group == "FallRx")
# 
# sa.piri2$predicted <- sa.piri2$predicted+0.01
# 
# sa.piri3 <- sa.piri1 %>% 
#   filter(group != "FallRx")
# 
# sa.piri3 <- rbind(sa.piri3, sa.piri2)

ggplot(sa.piri1, aes(x, predicted, color = group))+
  geom_line(linewidth = 1, show.legend = FALSE) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type")+
  theme(axis.text = element_text(size = 18), axis.title = element_text(size = 24))+
  labs(x = "Average leaf litter depth \n (standardized)",
       y = NULL)

ggsave(filename = "SA_PIRI_avgld.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)






sa.piri4 <- ggpredict(sa.m9s, terms = c("ba.ha_s", "Treat_Type"))
#again control and fall rx overlap - if I wanted to nudge/jitter

# sa.piri5 <- sa.piri4 %>% 
#   filter(group == "FallRx")
# 
# sa.piri5$predicted <- sa.piri5$predicted+0.01
# 
# sa.piri6 <- sa.piri4 %>% 
#   filter(group != "FallRx")
# 
# sa.piri6 <- rbind(sa.piri6, sa.piri5)


ggplot(sa.piri4, aes(x, predicted, color = group))+
  geom_line(linewidth = 1, show.legend = FALSE) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type")+
  theme(axis.text = element_text(size = 18), axis.title = element_text(size = 24))+
  labs(x = "Average basal area per hectare \n (standardized)",
       y = "Pitch pine stems \n (adjusted for time)")

ggsave(filename = "SA_PIRI_BA.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)





sa.piri5 <- ggpredict(sa.m9s, terms = c("other_s", "Treat_Type"))

ggplot(sa.piri5, aes(x, predicted, color = group))+
  geom_line(linewidth = 1, show.legend = FALSE) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type")+
  theme(axis.text = element_text(size = 18), axis.title = element_text(size = 24))+
  labs(x = "Counts of other saplings \n (standardized)",
       y = NULL)

ggsave(filename = "SA_PIRI_other.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)

Shrub Oak Models

Best Model AIC: 323.8

Start from the beginning knowing that both Region and Treat Type have 0 areas

tapply(sa.all4$Shrub_Oak, sa.all4$Region, summary)
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## $MA
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.000000 0.000000 0.006494 0.000000 1.000000 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   1.964   3.000   7.000 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2447  0.0000  6.0000
# no SO in LI,
tapply(sa.all4$Shrub_Oak, sa.all4$Treat_Type, summary) 
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.7094  0.0000  7.0000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.4902  0.0000  6.0000 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01786 0.00000 1.00000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0
# no SO in MowRX or SpringRx

rm(sa.all4)

looking at SO w/o spring rx or mowrx and w/ standardized variables

sa.all_noSMRX <- sa.all3 %>% 
  filter(Treat_Type != "SpringRx") %>% 
  filter(Treat_Type != "MowRx")

There is just not enough information in exposed mineral soil (so many 0s) for it to be useful. It messes up models, making them rank deficient. So just dropping that variable

so.sa1s <- glmmTMB(Shrub_Oak ~ Treat_Type + piri.ba_s + avgld_s + vegt_s + piri_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all_noSMRX,
                 ziformula = ~Region,
                 family = nbinom1)
# won't converge

#had to remove piri ba for model to converge
so.sa2s <- glmmTMB(Shrub_Oak ~ Treat_Type + avgld_s + vegt_s + piri_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all_noSMRX,
                 ziformula = ~Region,
                 family = nbinom1)
AIC (so.sa2s) #333
## [1] 333.03
#avgld
so.sa3s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + piri_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all_noSMRX,
                 ziformula = ~Region,
                 family = nbinom1)
AIC (so.sa3s) #331
## [1] 331.278
lrtest(so.sa2s, so.sa3s) #0.6
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + avgld_s + vegt_s + piri_s + other_s + 
##     ba.ha_s + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + piri_s + other_s + ba.ha_s + 
##     offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df Chisq Pr(>Chisq)
## 1  16 -150.51                    
## 2  15 -150.64 -1 0.248     0.6185
#piri
so.sa4s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all_noSMRX,
                 ziformula = ~Region,
                 family = nbinom1)
AIC (so.sa4s) #329
## [1] 329.3011
lrtest(so.sa3s, so.sa4s) #0.9
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + vegt_s + piri_s + other_s + ba.ha_s + 
##     offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + other_s + ba.ha_s + offset(l.TFT) + 
##     (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  15 -150.64                     
## 2  14 -150.65 -1 0.0231     0.8793
#other
so.sa5s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all_noSMRX,
                 ziformula = ~Region,
                 family = nbinom1)
AIC (so.sa5s) #327.3
## [1] 327.3061
lrtest(so.sa4s, so.sa5s) #0.9 - drop other
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + vegt_s + other_s + ba.ha_s + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + ba.ha_s + offset(l.TFT) + (1 | 
##     Site/Plot_No)
##   #Df  LogLik Df Chisq Pr(>Chisq)
## 1  14 -150.65                    
## 2  13 -150.65 -1 0.005     0.9436
#test BA
so.sa6s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all_noSMRX,
                 ziformula = ~Region,
                 family = nbinom1)
AIC (so.sa6s) #325.3
## [1] 325.3149
lrtest(so.sa6s, so.sa5s) #0.9
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + vegt_s + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + ba.ha_s + offset(l.TFT) + (1 | 
##     Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  12 -150.66                     
## 2  13 -150.65  1 0.0088     0.9253
# test veg
so.sa7s <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all_noSMRX,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.sa7s) #323.8
## [1] 323.7947
lrtest(so.sa7s, so.sa6s) #0.5
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  11 -150.90                     
## 2  12 -150.66  1 0.4798     0.4885
# add piri ba back in
so.sa8s <- glmmTMB(Shrub_Oak ~ Treat_Type + piri.ba_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all_noSMRX,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.sa8s) #325
## [1] 325.0044
lrtest(so.sa7s, so.sa8s) #0.4
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + piri.ba_s + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df LogLik Df  Chisq Pr(>Chisq)
## 1  11 -150.9                     
## 2  12 -150.5  1 0.7902      0.374
rm(so.sa1s,
   so.sa2s,
   so.sa3s,
   so.sa4s,
   so.sa5s,
   so.sa6s,
   so.sa8s)

Best model (no springrx or mowrx)

so.sa7s <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all_noSMRX,
                 ziformula = ~Region,
                 family = nbinom1)
summary(so.sa7s) #323.8
##  Family: nbinom1  ( log )
## Formula:          Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## Zero inflation:             ~Region
## Data: sa.all_noSMRX
## 
##      AIC      BIC   logLik deviance df.resid 
##    323.8    363.6   -150.9    301.8      264 
## 
## Random effects:
## 
## Conditional model:
##  Groups       Name        Variance       Std.Dev.  
##  Plot_No:Site (Intercept) 0.000000007557 0.00008693
##  Site         (Intercept) 0.109978958497 0.33163076
## Number of obs: 275, groups:  Plot_No:Site, 275; Site, 26
## 
## Dispersion parameter for nbinom1 family (): 0.773 
## 
## Conditional model:
##                   Estimate Std. Error z value            Pr(>|z|)    
## (Intercept)        -2.6952     0.2583 -10.435 <0.0000000000000002 ***
## Treat_TypeFallRx    0.8099     0.3689   2.195              0.0281 *  
## Treat_TypeHarvest   1.5834     1.2203   1.298              0.1944    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##                 Estimate   Std. Error z value Pr(>|z|)  
## (Intercept)       1.9034       0.7926   2.402   0.0163 *
## RegionLI         30.0378 1411573.3749   0.000   1.0000  
## RegionMA          2.0029       1.3083   1.531   0.1258  
## RegionME        -31.2661  873734.0264   0.000   1.0000  
## RegionNH         -1.0435       0.8939  -1.167   0.2431  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
so.sa7s_sr <- simulateResiduals(so.sa7s, n = 1000, plot = T) #passes

testResiduals(so.sa7s_sr)  #passes

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.032251, p-value = 0.9372
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9179, p-value = 0.902
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 275, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.005545455
## sample estimates:
## outlier frequency (expected: 0.00112727272727273 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.032251, p-value = 0.9372
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9179, p-value = 0.902
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 275, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.005545455
## sample estimates:
## outlier frequency (expected: 0.00112727272727273 ) 
##                                                  0

Pairwise

emmeans(so.sa7s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
##  Treat_Type response    SE  df asymp.LCL asymp.UCL
##  Control       0.722 0.187 Inf     0.435      1.20
##  FallRx        1.623 0.590 Inf     0.796      3.31
##  Harvest       3.518 4.244 Inf     0.331     37.42
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast          ratio    SE  df null z.ratio p.value
##  Control / FallRx  0.445 0.164 Inf    1  -2.195  0.0719
##  Control / Harvest 0.205 0.250 Inf    1  -1.298  0.3965
##  FallRx / Harvest  0.461 0.570 Inf    1  -0.626  0.8056
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale

none significantly different

Control vs FallRx close at 0.07

Now graph

Need to figure out better use of emmeans (potentially)

so.sa1.plot <- emmeans(so.sa7s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')

plot <- plot(so.sa1.plot, horizontal=FALSE, comparisons=TRUE, color = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15")) + theme_bw()

plot+
  theme_classic()+
  theme(panel.background = element_blank()) +
  theme(axis.text = element_text(size = 14), axis.title = element_text(size = 20), legend.position = "right", legend.text = element_text(size = 14), legend.title = element_text(size = 18)) +
  labs(x = "Shrub Oak Stems \n (corrected for time from treatment)",
       z = "Treatment Type",
       color = "Treatment Type")

More in harvest - does fire kill large scrub oak saplings? They are entirely missing from springrx and mowrx treatments

this is the same analysis, but with all treatment categories

There are no shrub oak sapling in springrx and mowrx treatments, so the above model is better - but just to have it.

so.all1 <- glmmTMB(Shrub_Oak ~ Treat_Type + piri_s + other_s + min_s + avgld_s+ piri.ba_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; function evaluation limit reached without convergence (9). See
## vignette('troubleshooting'), help('diagnose')
#won't converge

#drop min & piri ba

so.all2 <- glmmTMB(Shrub_Oak ~ Treat_Type + piri_s + other_s + avgld_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.all2) #335
## [1] 335.3884
so.all3 <- glmmTMB(Shrub_Oak ~ Treat_Type + other_s + avgld_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.all3) #333.5
## [1] 333.478
so.all4 <- glmmTMB(Shrub_Oak ~ Treat_Type + avgld_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.all4) #331.6
## [1] 331.5949
so.all5 <- glmmTMB(Shrub_Oak ~ Treat_Type + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.all5) #329.7
## [1] 329.725
so.all6 <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.all6) #327.8
## [1] 327.7947
so.all7 <- glmmTMB(Shrub_Oak ~ Treat_Type + piri.ba_s + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.all7) #329
## [1] 329.0044
rm(so.all1,
   so.all2,
   so.all3,
   so.all4,
   so.all5,
   so.all7)

so.all6_sr <- simulateResiduals(so.all6, n = 1000, plot = T) #passes

testResiduals(so.all6_sr) #passes

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.026629, p-value = 0.8718
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.88884, p-value = 0.844
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000522088353413655 ) 
##                                                   0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.026629, p-value = 0.8718
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.88884, p-value = 0.844
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000522088353413655 ) 
##                                                   0
emmeans(so.all6, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
##  Treat_Type response        SE  df asymp.LCL asymp.UCL
##  Control       0.410 0.1059479 Inf     0.247         1
##  FallRx        0.922 0.3350685 Inf     0.452         2
##  Harvest       1.998 2.4106105 Inf     0.188        21
##  MowRx         0.000 0.0000001 Inf     0.000       Inf
##  SpringRx      0.000 0.0000044 Inf     0.000       Inf
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                    ratio                   SE  df null z.ratio
##  Control / FallRx                0                    0 Inf    1  -2.196
##  Control / Harvest               0                    0 Inf    1  -1.298
##  Control / MowRx    10634534517938 16507512852342061056 Inf    1   0.000
##  Control / SpringRx    10155053986     1101249307522083 Inf    1   0.000
##  FallRx / Harvest                0                    1 Inf    1  -0.626
##  FallRx / MowRx     23902818311903 37103277047566393344 Inf    1   0.000
##  FallRx / SpringRx     22825109079     2475234065907171 Inf    1   0.000
##  Harvest / MowRx    51809117272110 80420978256787030016 Inf    1   0.000
##  Harvest / SpringRx    49473193394     5365044838234219 Inf    1   0.000
##  MowRx / SpringRx                0                 1486 Inf    1   0.000
##  p.value
##   0.1813
##   0.6927
##   1.0000
##   1.0000
##   0.9709
##   1.0000
##   1.0000
##   1.0000
##   1.0000
##   1.0000
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log scale